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The statistical property of acoustic emission (AE) events from a plunged granular bed is analyzed 
by means of actual time and natural time analyses. These temporal analysis methods allow us to 
investigate the details of AE events that follow a power-law distribution. In the actual time analysis, 
the calm time distribution and the decay of the event-occurrence density after the largest event (i.e., 
Omori-Utsu law) are measured. Although the former always shows a power-law form, the latter does 
not always obey a power law. Markovianity of the event-occurrence process is also verified using 
a scaling law by assuming that both of them exhibit power laws. We find that the effective shear 
strain rate is a key parameter to classify the emergence rate of power-law nature and Markovianity 
in the granular AE events. For the natural time analysis, the existence of self organized critical 
(SOC) states is revealed by calculating the variance of natural time Xk, where fcth natural time of 
N events is defined as Xfc = k/N. In addition, the energy difference distribution can be fitted by a 
g-Gaussian form, which is also consistent with the criticality of the system. 


I. INTRODUCTION 

Power-law distributions can be observed in various 
fields of natural and artificial phenomena [l|. Exam¬ 
ples include the distribution of incomes (Pareto’s law), 
the appearance frequency of English words (Zipf’s law), 
the number of meteorites impacting on a planet, and the 
number of species per genus in flowering plants 0. For 
seismic activity, one of the best known power-law rela¬ 
tions is Gutenberg-Richter (GR) law Q: 

G{Q) ^ Q-\ (1) 

where Q, G{Q), and 7 are the emitted energy per event, 
its occurrence frequency, and a characteristic exponent 
(positive constant), respectively. Recently, experiments 
and simulations of soft materials which are related to 
seismic activity have been performed. These studies have 
reported various power-law event-size distributions (e.g., 
sliding friction of gels Q and granular avalanches in sim¬ 
ulations i). Although this kind of power-law event-size 
distributions has also been discovered in many natural 
phenomena (e.g., forest fire areas floods [ 3 ], frag¬ 
ments Q, and Tsunami runup heights i), they are em¬ 
pirical laws and not fully understood in terms of their 
physical origin. 

The physical mechanisms determining the power-law 
exponent have long been studied. For instance, the ex¬ 
ponent of power-law size distribution in brittle fragmen¬ 
tation shows particular relations to the higher order mo¬ 
ment and system size [lol - [T3| . Besides, the exponent 
value depends on the dimensionality and the way of crack 
propagation Q. Also, stress relaxation mechanisms in 
plastic deformation can be classified by the power-law ex¬ 
ponent of stress drop distributions [l^. Although some 
fundamental aspects concerning the power-law size dis¬ 
tributions have been revealed as aforementioned, much 
deeper investigation associated with the appropriate clas¬ 
sification of the system is necessary for truly understand¬ 
ing the universality of the whole power-law distributions. 


In this study, acoustic emission (AE) burst events 
emitted from a plunged granular bed are particularly 
examined. In Ref. |I4| . the behavior of granular mat¬ 
ter has been studied using AE technique, in which the 
size distribution of AE burst events obeys power law. In 
the experiment, the power-law exponent (7 in Eq. ([T|)) 
varies depending on experimental conditions. Although 
the exponent value was related to the mode of defor¬ 
mation (brittle-like or plastic-like), the analysis has still 
been very qualitative in Ref. [l^. Thus, further detailed 
analyses are necessary to identify the underlying physi¬ 
cal mechanisms governing the power-law nature of gran¬ 
ular AE events. Such a deeper understanding of granular 
AE events might also provide a universal framework for 
various power-law distributions. Moreover, because the 
statistical behavior of dry granular matter is somewhat 
similar to that of seismicity [T^, the detailed study of 
granular behavior could be helpful to understand geo¬ 
physical phenomena as well. 

The most serious deficiency in the previous analyses of 
AE event-size distributions is a lack of temporal informa¬ 
tion. The power-law event-size distribution such as GR 
law usually neglects the time series of event occurrence. 
It only deals with the size of events whilst the events in¬ 
deed occur in the time series. In order to consider the 
temporal information, here we employ two analysis meth¬ 
ods: actual time and natural time analyses. 

In the actual time analysis, the amplitude of events is 
omitted in contrast to the analysis of event-size distribu¬ 
tions such as GR law. Then, the time interval between 
successive AE events called calm time (a.k.a. interoccur¬ 
rence time or waiting time) and the event-occurrence 
density are measured. The distribution of calm time 
has been found to be power law in many AE measure¬ 
ments (e.g., tensile failure experiment of paper sheets [T^ 
and volcanic rocks at Stromboli [13). For the event- 
occurrence density, the power-law decay of aftershock 
activity is known as Omori-Utsu (OU) law in seismol¬ 
ogy [1^, [13 ■ OU-like behavior is also observed in various 
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AE measurements (e.g., microfracturing in a compressed 
rock [13,[Si). 

If both these quantities (calm time and event- 
occurrence density) exhibit power-law distributions and 
each power-law exponent can be determined, we can eval¬ 
uate Markovianity of the event time series using these 
exponents. This powerful method to verify Markovianity 
was first developed for analyzing real seismic data, and 
non-Markov nature of the real seismic activity was re¬ 
vealed by Abe and Suzuki [H, [13| . In the current study, 
we are going to discuss the statistical property of gran¬ 
ular AE events by applying this method to the granular 
AE event data. 

The natural time analysis, on the other hand, discards 
the calm time information. It only uses the order of 
events and corresponding amplitudes. The idea of natu¬ 
ral time has also been proposed for the analysis of seis¬ 
micity [l^-H^. For an event time series comprising N 
events, natural time Xk serves as an index for the occur¬ 
rence of the fcth event and is defined as Xfc = k/N. Then 
the variance ki is computed as 

N / N \ 2 

Kl (2) 

k=l \k=l ) 


where pfc = Qkl^}2!i=\ Qi normalized released en¬ 

ergy Qk in the fcth event. The seismic electrical signals 
(SES) right before earthquakes tend to show a critical 
value Ki = 0.07 [2^. Similar ki behaviors can be con¬ 
firmed in some self organized critical (SOC) systems [IH- 
[S^, where the original concept of SOC was introduced 
by Bak et al. Furthermore, the AE signals from 

the deformed rock have exhibited the behavior similar to 
seismicity in terms of the natural time analysis [^ . 

In this study, the temporal properties of AE events 
emitted from a plunged granular bed are throughly in¬ 
vestigated through these analyses. The relation between 
the measured results and the previously obtained event- 
size distributions (GR law) [ij is also discussed on the 
basis of the experimental data. 


II. EXPERIMENT 

The experimental methodology and data used in this 
study are the same as those in Ref. [3. Glass beads 
(grain diameter d = 0.4, 0.8, or 2.0 mm) are poured into 
a cylindrical plexiglass container. A steel sphere (radius 
r = 5, 10, or 20 mm) is then penetrated into the granu¬ 
lar bed. The penetration speed is fixed as ?; = 0.5, I.O, 
or 5.0 mm s“^. The top surface of the granular bed is 
open to the atmosphere and any confining pressure is not 
applied to the bed. An AE sensor (NF AE-900s-WB) 
is buried and fixed in the granular bed to capture AE 
events created by the penetration. The AE sensor is a 
piezoelectric transducer which converts dynamic motions 
(e.g., ultrasonic elastic wave) into electric signals [s^ . 


Because the AE signals are very weak, they are ampli¬ 
fied by an amplifier (NF AE-99I3) and a discriminator 
(NF AE-9922). The sampling rate of the AE data is I 
MHz. Three experimental realizations for each setup of 
experimental conditions are performed to check the re¬ 
producibility. In general, granular behaviors have strong 
memory effect and history dependence [38l| . The sphere 
penetration must remain its memory in the granular bed 
such as the force chain structure in the tapped granu¬ 
lar bed [ 13 ]. Thus, the fresh granular bed is deposited 
before every experimental run to erase the memory of 
penetration. 

Let us summarize the experimental results briefly. A 
raw data example of AE signals A in volt as a function 
of time is shown in Fig. [Ija). While the origin of time in 
Fig. UDa) is defined by z = 0 (z is the penetration depth 
of the intruder) as similar to Ref. [l3 |. it is arbitrary in 
the current study. Actually, the origin of time t = 0 will 
be defined later by the main event. The corresponding 
experimental conditions are d = 0.8 mm, r = 10 mm, 
and V = 5.0 mm s“^. From now on, we use these experi¬ 
mental conditions for subsequent plots shown in this pa¬ 
per unless otherwise noted. Because all AE events show 
typical attenuating oscillation with short decay time as 
depicted in the inset of Fig. [Hb), each AE event can be 
picked up from the measured AE signals using a thresh¬ 
old value Ath = 0.06 V and deadtime tdead = 300 /rs [l3 |. 
The total number of AE events identified by this method 
ranges from 10^ to 10^. Figure [IJb) shows average power 
spectra of AE events for several experimental conditions. 
The dominant frequencies, 20 kHz and 75 kHz, seem to 
be independent of the experimental conditions and might 
result from complex coupling of both the AE device and 
the granular media. For more details of the experimental 
setting, see Ref. M- 


III. ACTUAL TIME ANALYSIS 


A. Calm time distribution 


First, we focus on the analysis of calm time t which 
corresponds to the time interval between two successive 
events. Specifically, the frequency distribution of calm 
time P{t) is measured. The occurrence time of each 
event is determined by the moment at which the signal 
amplitude exceeds the threshold value Ath- In Fig. [2Da), 
an example of P(t) distribution obeying the power-law 
form is presented. All other P{t) distributions also ex¬ 
hibit power-law forms. Namely, P(t) always obeys 


P{t) 


dN{T) 

Ntot 


t--(m-hi) 


( 3 ) 


where dNir) and Ntot are the number of events with 
calm time r intervals and the total number of events in 
the time series, respectively. The minimum (binning) 
timescale used in dfV(r) measurement is 1 ms. Since 
the number of events is limited and the finite size effect 
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must be considered carefully, here we employ the method 
of maximum likelihood estimation (MLE) |l|, Slj to 
determine the exponent value. MLE determines the ex¬ 
ponent to maximize the likelihood function L defined as 


L(1 -H/r|r) 


N 

n Ik _ 


(4) 


where Tk is the calm time of the fcth event normalized 
to the minimum value and C(1 + m) is the Riemann zeta 
function expressed as I2k=i To apply MLE, the 

exponent fi must be greater than zero. Otherwise, the 
Riemann zeta function diverges. Because MLE enables 
us to avoid lar ge b ias error in many frequency distribu¬ 
tions (e.g., (dol I42I. Id^lL MLE is a preferred method for 
accurately estimating the power-law exponent. The de¬ 
termined /i value in the data of Fig. ^a) is 0.82. To 




Frequency (kHz) 


FIG. 1. (Color online) (a) Example of the AE signals during 
the penetration of a sphere into a glass bead bed. (b) Aver¬ 
age power spectra of AE events produced per penetration for 
some experimental conditions. The inset shows an example 
of typical attenuating oscillation with short decay time. 


check the validity of the determined exponent from an¬ 
other aspect, the frequency distribution of r with the 
logarithmic bins using a constant rate is shown in 
the inset of Fig. UKa). Although the data in the inset 
are scattered, the global trend can be reproduced by the 
power-law fitting determined by MLE (p, = 0.82). 

For the obtained power-law exponent, fi depends on 
the ex per imental conditions. According to the previous 
study [3, the exponent of event-size distributions (GR 
law) 7 mainly depends on the grain size d. In contrast, 
the exponent of the calm time distribution P{t) is sen¬ 
sitive to various parameters. Figure [2](b) shows the v 
dependence of the exponent /i, where v and /r are posi¬ 
tively related. This means that the faster the penetration 
speed is, the more the relative frequency of shorter calm 
time increases. 


B. Event-occurrence density distribution 

Second, the number of events occurring per unit time 
after the main event S{t) is considered. S{t) is defined 

by 


S{t) = 


dN'jt) 

dt 


t-P, 


(5) 


where p is a characteristic exponent, t stands for the time 
elapsed after the main event, and dN'{t) is the number 
of events occurring in the short time interval between t 
and t + dt. We call P{t) the event-occurrence density. 
Because the time t = 0 should be defined by the moment 
of the main-event occurrence, it is necessary to locate 
the main event. In usual seismic activity, the main event 
(mainshock) can be located by subsequent smaller after¬ 
shocks [ij. Thus, in this work, we employ the largest AE 
event in the time series as a main event. Then, the event- 
occurrence density after the main event is measured from 
the experimental data. For real seismic activity, the uni¬ 
versal power law is well known as OU law, which states 
that S(t) decays with time following a power-law rela¬ 
tionship expressed on the right-hand side in Eq. ([5]). The 
exponent p usually ranges in 0.8 — 1.5 for real seismic ac¬ 
tivity [3^. 

However, the current experimental result shows various 
behaviors. Some dN'(t) distributions follow power-law- 
like decay (inset of Fig. EJa)), others seem to be almost 
constant (inset of Fig. [S^b)). The former and latter cor¬ 
respond to the nonsteady and steady processes, respec¬ 
tively. The nonsteady process means the existence of 
the main event followed by power-law-like decay of sub¬ 
sequent aftershock-like AE events. In other words, the 
meaningful main event cannot be identified in the lat¬ 
ter (steady) case. Actually, even in the power-law (non¬ 
steady) case, the scaled range is not very wide (only one 
order of magnitude) as can be seen in the inset of Fig. [S] 
This implies that it is not very easy to confirm a clear OU 
law (power-law decay of aftershocks) in the granular AE 
events. Although the decay range is limited, the inset of 
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FIG. 2. (Color online) (a) Calm time distribution dN{T) 
obeying a power law (Eq. (|3])). The inset shows dN{T) with 
logarithmically increasing bins, where the slope corresponds 
to —/i, not —(l+/i). The bin widths are given by 0.001(\/2)" s 
(e.g., 0.001, O.OOlv^, 0.002 ... s). Red broken lines in both 
plots show the power law with /r = 0.82 which is obtained 
by MLE method, (b) The exponent of the calm time distri¬ 
bution ^ as a function of the penetration speed v. The 
value depends on various experimental conditions in contrast 
to 7 in Eq. 0. The data of calm time used in this figure are 
obtained after the main event defined in Sec. ME 


Figs. m a) and (b) are significantly different. To discrim¬ 
inate these two phases and discuss the statistical prop¬ 
erty from the aspect of Markovianity, here we assume the 
power-law form for dN'(t). The estimated power-law ex¬ 
ponent values will play a crucial role to characterize the 
AE event statistics. 

In contrast to calm time distributions that show clear 
power-law behavior, care must be taken in precisely de¬ 
termining the fitted p value for event-occurrence den¬ 
sity distributions. Because most of fitted p values seem 
less than 1, the method of MLE cannot be used due 
to the divergence of Riemann zeta function. Instead, 
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FIG. 3. (Color online) Power-law histograms with logarith¬ 
mic bins of the event-occurrence density decay for (a) a non¬ 
steady (OU-like) AE event and (b) a steady (constant) AE 
event. Since the main plots show the logarithmically binned 
data, OU law is expressed as dN'{t) = t^~’^dlnt. Here, the 
range dt = 0.1('\/2)° — 0.1('\/2)® s is used. The least-squares 
method is used for determining the p value. The insets show 
the power-law histograms dN'{t) ~ S{t)dt with the constant 
linear binning dt = 0.3 s 


logarithmically increasing bins are employed here. Un¬ 
equal bin widths are usually used to obtain a more ho¬ 
mogeneous number of data per bin than that with a 
constant bin width, which can reduce statistical errors 
in t he p ower-law tail due to the poor number of sam¬ 
ples [43,|43- Here, the bin widths are given by 0.1(-\/2)" s 
(e.g., 0.1, 0.1-\/2, 0.2 ... s). FiguresISKa) and (b) show the 
histogram of the number of aftershock-like events dN'{t) 
after the main event obtained by using the logarithmic 
bins. By applying the logarithmic binning, the power-law 
exponent of the data plot varies as. 


dN'jt) 

dint 


s{t)t - t^-p. 


( 6 ) 


Namely, the slopes in the main plots of Fig. [3] corre¬ 
spond to I — p. Although the data considerably scatter 
in Fig. |3Ka), we assume the power-law behavior at least 
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in the range of 0.1 < t < 3 s. Then, the Markovianity of 
the event series after the main event can be evaluated as 
discussed later in Sec. lIII Cl Using this procedure, all the 
data are fitted by power-law forms. The fitted p values 
for Figs. a) and (b) are 0.39 and —0.02, respectively. 
The corresponding slopes are also shown in the insets as 
well. 

In Fig. m the histogram of the fitted p values for all 
the dataset is presented. As can be seen in Fig. 01 the 
histogram shows a bimodal distribution and the valley 
around p = 0.3 seems to discriminate two phases. Thus, 
we define the case of p > 0.3 as nonsteady (OU-like) 
states. One can also confirm a small portion of popula¬ 
tion in the negative p regime. This might be interesting 
behavior, possibly indicating the precursor of the large 
event. However, here we only focus on the power-law 
decay corresponding to the OU-like state. 

As discussed thus far, some experiments show the OU- 
like behavior and others do not. What is the most im¬ 
portant parameter determining the behavior of the event- 
occurrence density? To answer this question, we study 
some parameter dependencies and find that the effective 
shear strain rate v/r (the penetration speed divided by 
the radius of a penetrator) is relevant to character¬ 
ize the state. Specifically, the emergence rate of OU- 
like states is shown as a function of v/r in Fig. [SJ The 
emergence rate of OU law is defined by the ratio of the 
experimental realizations showing p > 0.3 to the total 
experimental realizations with the identical v/r. Since 
the number of experimental realizations with the identi¬ 
cal v/r is not constant, we employ the normalized occur¬ 
rence ratio to estimate the emergence frequency. One can 
confirm that the emergence rate of OU law (nonsteady 
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FIG. 4. (Color online) Histogram of the fitted p value showing 
a bimodal structure. To distinguish the nonsteady (OU-like) 
distributions from the steady distributions, the value at the 
valley between two peaks (p = 0.3) is used (vertical broken 
red line). 


state) abruptly increases in v/r > 0.25 s“^. Although 
the complete reproducibility of OU-like behavior is not 
established even in the range of v/r > 0.25 s“^. Fig. [5] 
implies that v/r = 0.25 s“^ is the marginal value between 
the steady and nonsteady (OU-like) regime. 

Next, to verify the v/r dependence of p, the averaged p 
value of OU law is plotted as a function oiv/r in Fig. El 
In the OU-like regime {v/r > 0.25 s“^), the p values 
show a roughly constant value, p ~ 0.45, in contrast to 
the steady regime suggesting p ~ 0.3, which might result 
from the defined marginal value between steady regime 
and OU-like regime (p = 0.3). Since the emergence rate 
of OU-like behavior is very low in v/r < 0.25 s“^, the 
value is sensitive to the threshold. However, the data 
errors in Fig. El are considerably large, thus the difference 
between p = 0.3 and p = 0.45 is not very clear. 


C. Markov scaling 

Using two exponents p and p dehned in Eqs. m and 
(El) , Markovianity of the event time series after the main 
event can be discussed in terms of a scaling law. The 
specific values used in the evaluation of the Markov scal¬ 
ing are p (Fig. 0]) and p (Fig. El) computed from the data 
after the main event. The scaling law used here was 
originally developed and applied to real seismicity. As a 
result, non-Markov nature of earthquake aftershocks was 
reported [13, [1^. Here, we slightly expand this scaling 
law. If a process of event occurrence is Markovian, the 
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FIG. 5. (Color online) Effective shear strain rate (v/r) depen¬ 
dence of the emergence rate of OU-like behavior. The larger 
shear strain rate {v/r > 0.25 s“^) tends to cause more OU-like 
behaviors than the small shear strain rate {v/r < 0.25 s“^). 
We express these two distinctive regimes as steady and OU- 
like {nonsteady) regimes. 
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following equation holds (46|: 


the resultant Laplace transformation becomes 


S(t)=P{t)+ [ P{t-t')S{t')dt', (7) 

JQ 

which can be derived from the Kolmogorov forward equa¬ 
tion [d^. Then, the Laplace transformation of Eq. © 
yields 


£[P]{s)^l-/3s, (14) 

for a small limit of s, where /3 is a positive constant. 
Using Eqs. m and m, Eq. ([8]) results in 

p = 0. (15) 


C[S]{s) 


l-£[P](s)’ 


( 8 ) 


where £[/](s) = f{t)dt. Here, we assume that 

both P{t) and S{t) decay following the power-law forms 
as written in Eqs. © and © for a large value of t. Then, 
the Laplace transformations of P{t) and S{t) result in 
different expressions depending on the ranges of the ex¬ 
ponents. If the exponents fj, and p are in the ranges 


0 < p < 1 and 0 < p < 1, (9) 

the Laplace transformations of P{t) and S{t) behave as 
£[P](s) - 1 - (10) 


( 11 ) 

s ^ 

for a small limit of s, where a is a positive constant. 
Substituting Eqs. TO and (|lll) into Eq. ®, we obtain a 
simple scaling relation : 

p + p = l- ( 12 ) 

However, if the exponent p lies in the range 

1 < p < 2, (13) 


Note that this result is inconsistent with the assumed 
range of Eq. this implies that Eq. m is only an 
asymptotic solution. In the case of /i > 2, the scaling 
relation fulfilled in a Markov process does not exist. 

Equation (TT^ indicates the criterion for the Marko- 
vianity in the OU-like (nonsteady) event time series in 
the range of Eq. ®. To verify Markovianity of the gran¬ 
ular aftershock-like AE events showing OU-like behavior, 
p -|- /r as a function oi v/r is depicted in Fig. 0 The bro¬ 
ken red line in Fig. [7] indicates the Markov scaling law 
(Eq. (fT^ L One can confirm that some data (triangles) 
show clearly large values {p + p > 1 ) that are somewhat 
similar to earthquake aftershocks [22, l23| ■ In the current 
experimental result, however, the other data points (cir¬ 
cles) distribute around p-|-p~ 1 (orp-|-p< 1). This 
means that there might be a certain parameter range in 
which the event time series can be regarded as a Markov 
process. This result is contrastive to the real seismic ac¬ 
tivity that always shows p + p > 1, i.e., a non-Markov 
property. 

Actually, the effective shear strain rate v/r is not a 
unique parameter to characterize the system. To de¬ 
scribe the global behavior of all granular AE events in 
this study, the normalized grain size d/r is also used here. 
Indeed, d/r can be a characteristic parameter to classify 
the event-size distribution (GR law) [T^. d/r = 0.04 
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FIG. 6. (Color online) Effective shear strain rate {v/r) de¬ 
pendence of the fitted p value. The broken red line represents 
the level of p = 0.45 while the solid black line refers to the 
marginal value p = 0.3. The p value obtained in OU-like 
regime {v/r > 0.25 s“^) is almost independent of v/r. 


FIG. 7. (Color online) p -I- p as a function of the effective 
shear strain rate v/r. The broken red line indicates the level 
p -I- p = 1 (Markov scaling law). Circles distribute around 
p -I- p ~ 1 (or p -|- p < 1) while triangles clearly beyond the 
value. 























7 


is regarded as a marginal value between brittle-like and 
plastic-like behavior. As discussed in the last subsec¬ 
tion, on the other hand, vjr = 0.25 s“^ is considered 
as a marginal value between the steady and nonsteady 
(OU-like) states. Therefore, we make a phase diagram as 
shown in Fig. [ 8 l where the vertical and horizontal axes 
indicate the effective shear strain rate vjr and normal¬ 
ized grain size d/r, respectively. Although it is difficult 
to draw clear phase boundaries, there might be several 
classes of behavior in the phase diagram. In the large vjr 
and large djr regime, the nonsteady (OU-like) event time 
series, which is shown by triangles, can be frequently ob¬ 
served. In this OU-like regime, the value of |p -I- ^ — 1| is 
indicated by the gray scale in each triangles as long as 
is less than I (i.e., the range defined by Eq. ®). Namely, 
the dark triangle means that the process is non-Markov. 
The steady {p ~ 0) cases are represented by circles. The 
open and filled circles correspond to the steady Markov 
(1 < /i < 2, p ~ 0) and steady non-Markov (0 < /i < I 
or 2 < /i, p ~ 0 ), respectively. 
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FIG. 8. (Color online) Phase diagram of the temporal statis¬ 
tics of granular AE events, x- and y-axes represent d/r and 
v/r, respectively. The meaning of each symbol is as follows. 
The triangles represent the cases in which nonsteady (OU- 
like) behaviors can be observed. Note that triangles does 
not mean the complete reproducibility of OU-like behavior. 
As mentioned in the text, three experimental realizations for 
each experimental conditions were carried out. When two of 
three realizations show OU-like behavior, the triangle mark 
is used. The gray scale in triangles indicates the value of 
\p + fi — 1\ as shown in the legend. The open circle indicates 
the steady Markov process (l</i<2,p~0), while the filled 
circles indicate the steady non-Markov process (0 < p < 1 or 
2 < p, p ~ 0). The vertical line d/r = 0.04 is the boundary 
between brittle-like and plastic-like regimes [ 13 ] . The hor¬ 
izontal dashed line v/r = 0.25 s“^ is considered to be the 
characteristic shear strain rate above which the OU-like be¬ 
havior can often be observed. 


IV. NATURAL TIME ANALYSIS 


A. Variance ki 


Next, the natural time analysis is applied to the iden¬ 
tical dataset. In the definition of natural time, calm time 
is completely neglected. Instead, the order and magni¬ 
tude of events are utilized to characterize the event statis¬ 
tics. By simply applying Eq. m to the dataset, we can 
readily calculate ki. Here we use the squared maximum 
amplitude of each AE event for its released energy value 
Q k H- Some examples of Ki{k) are shown in Fig. |9l 
As expected, ki looks fluctuating around ki ~ 0.07 in 
some data, which indicates the criticality of the system. 
However, this tendency is not very universal. For in¬ 
stance, the asymptotic value of Ki{k) in Fig. EKc) seems 
to be different from 0.07. To verify the ki behavior in 
more detail, the probability density function (PDF) of ki 
is computed by the following procedure [H, [s^. First, 
the Ki value is computed from the natural time windows 
for 6 — 40 consecutive events. Second, this process is 
performed for all events by scanning the whole data set. 
Figure [To] shows the PDF computed from the data used 
in Fig. 0 The most probable value Ki^p estimated by 
the peak location of the PDF depends on experimental 
conditions. Particularly, the grain diameter d seems to 
be an important parameter. This d-dependent tendency 
is similar to the behavior of 7 in Eq. ([T|) (event-size dis¬ 
tribution exponent) M- 
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FIG. 9. (Color online) Evolution of ki as a function of the 
number of events k for various experimental conditions of (a) 
d = 2.0, (b) 0.8, and (c) 0.4 mm from the top to the bottom, 
respectively. The other parameters are fixed as r = 10 mm 
and V = 5.0 mm s“^. The evolutions of each ki are shown 
up to fc = 3000. The horizontal broken black lines represent 
Ki = 0.07 indicating the criticality of the system. The top (a) 
and middle (b) data look fluctuating around 0.07 while the 
data in (c) is clearly offset. 





















B. Returns distribution 

Although the measurement of ki is easy and useful to 
briefly characterize the critical state in the time series, 
it is in general not sufficient to evaluate the criticality of 
the event time series. Caruso et al. [3 have suggested 
another way for characterizing the criticality in the sys¬ 
tem. 

The variable returns x{k) are defined as x{k) = Qk+i — 
Qk- Namely, x{k) corresponds to the energy difference 
between two successive events. Here the returns are 
normalized to the mean {x) and standard deviation u, 
5 = {x—{x))/a. Then the PDF of ^ is calculated for each 
experimental condition (e.g., the inset of Fig. ni). As a 
result, it is clarified that the functionaf form of the PDF 
is independent of experimental conditions. Therefore, to 
obtain better statistics, the whole data of the normal¬ 
ized returns under various experimental conditions are 
merged into a single PDF. The entire PDF as a function 
of 5 is shown in the main plot of Fig.[TTl in which the PDF 
has much broader tails than normal Gaussian distribu¬ 
tion /(^) = exp(—^^/2)/-\/2^ (broken black curve). In¬ 
stead, q-Gaussian form/(^) = — — 

can fit the data well (red curve), where Aq and Bq are 
constants [i^. The computed q value is q = 1.9. If the 
data of a single experimental run are used for the PDF 
as shown in the inset of Fig. 1111 the agreement between 
q-Gaussian and data is limited particularly in the tail 
part. However, the qualitative tendency of the distribu¬ 
tion is basically identical. Because the g-Gaussian PDF 
of ^ can be observed in the SOC model satisfying both 
the power-law event-size distribution and the finite size 
scaling [d^, this result provides the supportive evidence 



FIG. 10. (Color online) PDF of the ki value in the exper¬ 
imental conditions same as Fig. The color code is iden¬ 
tical to that in Fig. [O] The vertical broken line indicates 
= 0.07. The peak values of the blue {d = 2.0 mm) and red 
{d = 0.8 mm) curves are around 0.07. 


for the presence of the critical state. 


C. Relation between ki and 7 

Finally, we consider the relation between natural time 
analyses and event-size distributions. As shown in Figs. [9] 
and m the Ki value is mainly affected by the grain di¬ 
ameter d. Additionally, the exponent of power-law event- 
size distributions 7 (Eq. [T]) also depends on dM- Thus, 
Ki and 7 might show a certain relation. To check the 
relationship, ki vs. 7 is plotted in Fig. [T^l As expected, 
a correlation between them can be confirmed. The sym¬ 
bols (or colors) in Fig. [T^] indicate the difference of de¬ 
formation mode: brittle-like mode characterized by the 
smaller 7 (blue circles) and plastic-like mode character¬ 
ized by larger 7 (red triangles). The solid black curve 
in Fig. [ 12 ] provides a fitting by exponentially asymptotic 
function Ki^p = 0.083 — The asymptotic value 

coincidentally agrees with Ku = 0.083 (~ 1/12) of a uni¬ 
form distvihution [25l - [^ [^. In Fig.[T^ Ki = 0.07 seems 
to be satisfied in the experimental conditions of brittle- 
like behavior. Put simply, the brittle-like regime shows 
Ki ~ 0.07 while ki approaches 0.083 by increasing 7 , i.e., 
in more plastic-like (flowing) regime. The condition that 
Ki approaches 0.083 means that the critical state is not 
established in the time series of AE events. 





FIG. 11. (Golor online) PDF of the normalized returns ^ 
computed from the whole AE events in all experimental con¬ 
ditions. The data obtained in the current experiment is bet¬ 
ter fitted by g-Gaussian (solid red curve) than normal Gaus¬ 
sian (broken black curve). The inset shows an example of 
single PDF in the experimental conditions of d = 0.8 mm, 
r = 10 mm, and v = 5.0 mm s“^. 
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V. DISCUSSION 

Thus far, various temporal analysis methods have been 
applied to the granular AE event data. Here let us discuss 
their physical meaning and relations. 

In the actual time analysis, power-law exponents for 
the calm time distribution and the decay of the event- 
occurrence density were measured. All of the experi¬ 
mental data show power-law forms for the calm time 
distribution. However, the obtained exponent value /i 
significantly depends on the penetration speed v. This 
tendency is in contrast to the case of power-law event- 
size distributions (GR law), in which the exponent 7 is 
mainly determined by the grain size d and almost in¬ 
dependent of the penetration speed v. For the event- 
occurrence density, power-law behavior is not universal. 
The emergence rate of the power-law decay (OU-like be¬ 
havior) can be characterized by the effective shear strain 
rate v/r. Since the dimension of v/r corresponds to the 
inverse of time, this quantity indirectly represents a char¬ 
acteristic timescale in the penetration system. Addition¬ 
ally, the time range of OU-like behavior is equal to or less 
than the order of 10° s (e.g.. Fig. [3]), suggesting the relax¬ 
ation process occurs within the characteristic timescale 
(0.25“^ = 4 s). Therefore, it is natural that the tem¬ 
poral property such as the event-occurrence density can 
be sorted by v/r. In the range of the OU-like regime 
in the granular AE events {v/r > 0.25 s“^), the power- 
law exponent p shows an almost constant value p ~ 0.45 



FIG. 12. (Color online) The most probable value Ki,p lo¬ 
cated by the peak of the PDF as a function of the power-law 
event-size distribution (GR law) exponent 7 . The colors and 
symbols represent the difference between brittle-like regime 
(blue circles) and plastic-like regime (red triangles). The blue 
circles tend to distribute around 0.07. The asymptotic value 
estimated by the fitting is approximately 0.083, which coin¬ 
cidentally corresponds to ki of a uniform distribution [H- 

[13, m. 


independent of v/r (Fig. [ 6 ]). This p value is less than 
the typical value for real seismicity and AE events from 
the microfracturing of rocks, p ~ 1 M- The reason for 
this discrepancy remains unsolved. Perhaps, this differ¬ 
ence originates from the peculiar nature of deformation 
in bulk granular matter. 

For the natural time analysis, on the other hand, the 
variance ki is related to the exponent of power-law event- 
size distributions (GR law) 7 . Since the exponent 7 is 
mainly determined by the grain size d dl , the Ki in gran¬ 
ular AE event time series is also related to the grain size 
d. Additionally, in the natural time analysis, the interval 
time between events (calm time) is completely neglected 
and only the order and amplitude of events are used. This 
is the reason why the geometrical parameter such as d or 
d/r becomes more essential than the temporal parameter 
u or u/r in the natural-time-related analysis. These two 
analysis methods (actual time and natural time analyses) 
are complementary to each other. 

As a matter of fact, the correlation between Ki^p and 
7 similar to Fig. IT^ has also been observed in the arti¬ 
ficially randomized (shuffled [ 13 ) event series data [30| . 
However, the agreement remains qualitative; the specific 
value ranges are slightly different between the random¬ 
ized event series and the current experimental data. This 
difference might come from the effect of memory among 
events as the randomized data do not have memory. 

The statistical behavior of granular AE events be¬ 
comes similar to that of real seismicity when both d/r 
and v/r are in the appropriate ranges: d/r > 0.04 and 
0.25 < v/r < 1.0 s“^ (Figs. [7l El and [12] and Ref. H)- 
In this regime, the values of 7 , ki, and p-\- p show similar 
values between the granular AE events and real seismic 
activity. Note that, however, the specific values of p and 
pL are different between the granular AE events and real 
seismicity. Furthermore, this coincidence in character¬ 
istic quantities might not directly mean the correspon¬ 
dence of underlying physical mechanisms. 

For instance, OU law for real seismic activity can be 
considered as the relaxation process after the mainshock. 
When the penetration speed is rapid, the available time 
for the relaxation becomes relatively short (might be in¬ 
sufficient) in general. Thus, it is difficult to see the re¬ 
laxation in the large v/r regime. On the other hand, 
the large v/r is necessary to reproduce OU-like behavior 
in the granular AE experiment. This implies the large 
v/r might result in the large shear field which causes 
a number of aftershock-like AE events. These two ef¬ 
fects will compete with each other and the qualitative 
tendency seems to be opposite. As a result, the com¬ 
plex statistical properties are observed. To simplify the 
problem, the controlled experiment in a different geomet¬ 
rical setup should be performed. The current experimen¬ 
tal setup actually originates from the previous researches 
concernir^ the slow-penetration drag force in granular 
matter [ 3 , [ 3 ) 113 ■ The simpler setup such as a simple 
shear, which can mimic the fault slip, might be better for 
future studies. 
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VI. CONCLUSION 

The statistical property of granular AE events was in¬ 
vestigated using two approaches: actual time analysis 
(Sec. mil) and natural time analysis (Sec. ITVT) . In the 
actual time analysis, the calm time distribution always 
shows a power-law form while OU-like behavior can only 
be frequently observed in the range of vjr > 0.25 s“^. 
In the OU-like (nonsteady) regime, non-Markov behav¬ 
ior is observed in a particular vjr regime. However, the 
steady (not OU-like) behavior is actually dominant in 
the granular AE events. To control the emergence of OU- 
like behavior and Markovianity, the appropriate tuning of 
vjr (inverse of the timescale) is necessary. Natural time 
analyses revealed that the ki value distributes around 
0.07 — 0.083. In addition, q-Gaussian fits the distribu¬ 
tion of the returns (i.e., energy difference among events). 
This result supports the SOC state of the system. The ki 
value can be related to the exponent of power-law event- 
size distributions (GR law) 7 . This means that the grain 
size d is the important parameter for ki because 7 is di¬ 


rectly related to d [l3|. We find that m ~ 0.07 can be 
established in the brittle-like regime {d/r > 0.04). 

In summary, statistical properties of seismic activity 
can be mimicked by the granular AE events in the range 
d/r > 0.04 and 0.25 < v/r < 1.0 s“^. Although the 
current experimental system is different from the mi¬ 
crofracturing of rocks and geological scale phenomena, 
the AE data obtained from a plunged granular matter 
exhibits some similarities with geological scale phenom¬ 
ena like earthquakes in terms of actual time and natural 
time analyses. 
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